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A METHOD  FOR  THE  NUMERICAL  SOLUTION  OF  A 
PARTICULAR  SET  OF  COUPLED  ORDINARY 
NONLINEAR  DIFFERENTIAL  EQUATIONS 


INTRODUCTION 


The  mathematical  modeling  of  certain  physical  phenom- 
ena can  be  treated  by  many  basic  methods.  Certain  cases 
may  lead  to  systems  of  nonlinear  differential  equations 
which  can  only  be  investigated  by  approximate  methods.  This 
paper  discusses  a system  of  equations  of  the  form 

AQ  = BG  - CQ2  - DQ  - EQ, 

where  A^/CjD^and  E are  square  coefficient  matrices  and 

Here  is  the  column 


A method  by  which  numerical  solutions  may  be  effected  is 
also  discussed. 


G,Q.Q,Q  and  are  column  vectors, 
vector 
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PARTICULAR  SET  OF  EQUATIONS 


Consider  the  system  of  differential  equations  written 
in  matrix  form  as  below, 

AQ  = BG  - CQ2  - DQ  - EQ,  (1 

where  the  matrices  and  vectors  are  real  and  A is 
nonsingular . 

If  A is  a square  nonsingular  matrix  then  it  possesses 
an  inverse,  A“l.  Multiplying  through  by  A-^  gives 

Q = A -iBG  - A_1CQ2  - A-1DQ  - EQ  (2a 

or 

Q = JG  - KQ2  - LQ  - MQ,  (2b 

where 

J = A-1B 
K = A-1C 
L = A_1D 
M = E. 

If  Q is  a column  vector  of  order  n,  then  Q can  be 
expressed  as 


ft 
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and  its  time  derivatives  as 


and 


(4) 


(5) 


If  A,  B,  C,  D,  and  E are  of  (n  x n)  order,  the  matrix 
multiplication  in  (2a)  will  yield  n ordinary  coupled  differ- 
ential equations  which  are  nonlinear  due  to  the  q?  terms  as 
shown  below, 


• • r- 

qi  = fl 

(t, 

...  qn, 
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n 
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• * * 
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.(6c) 

At  the  present  time,  a closed-form  solution  to  these 
equations  is  not  obtainable.  It  should  be  noted  that  for 
some  physical  problems,  the  coefficient  matrices  may  possess 
time-dependent  elements;  if  so,  equations  (6)  become  a set 
of  nonautonomous , nonlinear,  differential  equations. 
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. 
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A BRIEF  BACKGROUND  ON  THE  RUNGE-KUTTA  METHOD 


The  solution  of  nonlinear  differential  equations  in 
closed  analytic  form  is  not  usually  possible.  The  demands 
of  applied  science  make  it  necessary  to  obtain  some  insight 
into  the  nature  of  solutions  subject  to  given  boundary  con- 
ditions. Graphical  and  numerical  representation  of  the 
functions  are  the  usual  methods  employed. 1 

The  numerical  solution  of  differential  equations  by  a 
method  of  finite  integration  in  one  of  its  several  forms  is 
a favorite  tool  when  performed  by  means  of  digital  comput- 
ers. The  Runge-Kutta  method  is  used  perhaps  more  than  any 
other,  and  will  be  reviewed  here  briefly. 2 

Consider  a system  of  first-order  simultaneous  differ- 
ential equations.  For  simplicity,  choose  two  as  given  below 
(note  that  the  method  can  be  extended  to  n equations) . 


f(t,  y. 

z) 

(7) 

dz 

dt 

g(t,  y. 

z)  . 

The  estimate  for  the  function  y and  z at  some  later 
time  (t  + h)  can  be  computed  by  the  relations 

y(t)  = yQ  + Ay  = yQ  + 1/6  (A' z + 2A"y  + 2A"'y  + A4y) 

and  (8) 

z(t)  = zQ  + Az  = zq  + 1/6  (A' z + 2A'z  + 2A"'z  + A4z)  . 

The  increment  values  (A1)  are  computed  from: 

A'y  = f(tQ,  yQ,  zq)  h (9) 

A'  z = g ( t , y , z ) h 
o o o 

A"y  = f(tQ  + l/2h , yQ  + 1/2  A'y,  zQ  + 1/2  A'z)  h (10) 

A"z  = g(tQ  + l/2h,  yQ  + 1/2  A'y,  zq  + 1/2  A'z)  h 
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A"  'y  = f(tQ  + l/2h,  yQ  + 1/2  A"yo,  zq  + 1/2  A”z)  h 
A"'Z  = g(tQ  + l/2h,  yQ  + 1/2  A"y , zQ  + 1/2  A"z)  h 


(ID 


A4y  = f(tQ  + h,  yo  + A"'y,  zq  + A"  ’z)  h 
A4z  = g(tQ  + h,  yQ  + A"  ' y , zq  + A" ’z)  h. 


(12) 


If  the  increments  of  y and  z are  computed  in  the  order 
given , only  the  previously  computed  increments  are  needed  in 
each  step  in  the  computation  (i.e.,  the  increments  cannot  be 
computed  simultaneously) . 

This  method  can  be  applied  to  a system  of  higher-order 
equations  without  essential  change.  Consider  the  system  of 
second-order  equations 


and 


^ - f<t.  Y.  z,  ) 

dt2  dt  dt 


= g(t,  y,  z,  ) . 

dt2  dt  dt 


(13) 


Introducing  the  new  variables 

u = ^ and  v = dz 


dt 


dt 


(14) 


the  system  of  equations  given  by  (13)  becomes 


and 


— = f (t,  y,  z,  u,  v) 
dt 


dv  . . . 

— = g (tf  y,  z,  u,  v). 
dt 


(15) 


This  system  of  equations  can  now  be  numerically  inte- 
grated, using  the  method  described  earlier. 

Consider  the  system  of  equations  given  in  (6)  earlier. 
The  Runge-Kutta  integration  method  requires  that  the 
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highest-order  derivative  of  only  one  dependent  variable  be 
on  the  left  side  of  the  equality  in  each  equation,  since  the 
estimate  of  each  function  (see  equations  (9)  through  (12)) 
is  computed  in  a stepwise  fashion,  using  the  increment 
values  A1.  An  attempt  to  numerically  integrate  a system  of 
equations  with  more  than  one  highest-order  derivative  on 
either  side  of  the  equality  will  make  the  Runge-Kutta 
algorithm  unstable.  The  set  of  enuations  given  in  (6)  are 
now  of  a form  which  is  acceptable  to  the  Runge-Ku*-ta 
algorithm. 

Note  that  the  elements  of  the  coefficient  matrices  may 
be  time-dependent.  If  so,  they  can  be  computed  at  each  time 
increment  before  the  matrix  operations  take  place.  Most 
digital  computers  are  equipped  with  built-in  subroutines 
for  the  evaluation  of  matrices  and  determinants.  These  can 
be  incorporated  within  the  numerical  integration  loop  as 
the  integration  is  carried  on. 


EXAMPLE 

Consider  the  equation  of  motion  of  a tow.-d  body  movinq 
through  a nonaccelerating  viscous  fluid. * • 4 In  body 
coordinates,  the  equations  become,  in  mat nc  form, 

MSQ  - MSG  + H + T = 0 , (16) 

where 

Ms  = the  body  structural  mass  matrix 
Q = the  displacement  vector  of  the  body 
G = the  gravitational  vector 
H = the  hydrodynamic  force  vector 
T = the  towline  tension  vector. 

The  hydrodynamic  force  vector  (from  boundary  layer 
theory 5)  can  be  written  as 

H = MhQ  + D,  (17) 

where 

Mft  = the  body  hydrodynamic  mass  matrix 
= mhi j (see  appendix) 

D = the  viscous  force  vector. 
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Substituting  (17)  into  (16)  gives 

(Ms  + Mh)  Q = MG  - D - T.  (18) 

If  the  viscous  force  vector  is  considered  to  be  given  in  the 
form 

D = CxQ2  + C2Q,  (19) 

where  Ci  and  C2  are  viscous  coefficient  (i.e.,  lift,  drag, 
linear-damping  coefficients  or  stability  derivatives) 
matrices  for  the  body  in  question,  substitution  into  (12) 
gives 

(Ms  + Mh)  Q = MSG  - C i(i2  - C2Q  - T,  (19) 

which  is  of  the  same  form  as  (1)  discussed  previously. 

For  a body  with  six  degrees  of  freedom,  there  are 
6 x 2 = 12  equations  which  describe  the  body  motion  (i.e., 
acceleration  and  velocity)  and  the  necessary  auxiliary  re- 
lations which  must  be  found  to  compute  the  towline  tension. - • 4 

Using  the  method  described  in  the  first  section,  these 
equations  can  be  numerically  integrated  when  put  into  the 
form  as  given  by  (6).  Note  that  the  elements  of  the  hydro- 
dynamic  mass  matrix  may  be  time-  or  amplitude-dependent . 6 , 7 
For  this  case,  they  can  be  recomputed  at  the  end  of  each 
time  increment  in  the  numerical  integration  before  the 
matr..x  operations  in  (2)  are  performed. 
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SUMMARY 


It  has  been  noted  that  a particular  set  of  nonlinear 
differential  equations  may  be  reduced  to  a form  such  that 
numerical  integration  of  the  equations  may  be  effected  using 
a Runge-Kutta  method.  One  example  of  a physical  system  to 
which  a set  of  equations  of  this  form  may  apply  has  been 
provided  to  illustrate  the  problem.  It  is  the  intent  of 
this  papef  to  point  out  one  approximate  method  of  solution 
to  these  equations  as  an  aid  to  anyone  confronted  by  a 
similar  problem. 
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APPENDIX 

HYDRODYNAMIC  MASS 


For  the  case  of  a body  moving  through  a fluid,  the 
effects  on  body  motions  due  to  hydrodynamic  inertia  (or 
hydrodynamic  mass)  terms  must  sometimes  be  accounted  for  by 
the  mathematical  model  used.  For  a body  moving  in  an  un- 
steady manner,  the  fluid  disturbance  due  to  the  body  motion 
extends  in  decreasing  amplitude  to  infinity.  From  a system 
of  reference  moving  with  the  body,  the  effect  can  be  con- 
sidered as  localized  to  a finite  volume  of  fluid  surrounding 
the  body.  The  body  acts  as  if  an  added  mass  of  fluid  is 
moving  with  it.  If  the  body  is  subjected  to  an  accelera- 
tion, not  only  must  the  body  mass  be  accelerated  but  also 
the  added  mass  of  fluid. 8 

Writing  Newton's  Second  Law  for  the  body  gives 

F = (m  + m,  ) a . 
h 

Here  m is  the  mass  of  the  body  and  m is  the  added  mass  or 
hydrodynamic  mass  of  the  body. 

The  hydrodynamic  mass  is  expected  to  be  proportional 
to  the  fluid  density,  body  size  (i.e. , its  volume),  and  its 
orientation  while  in  motion  through  the  fluid.  For  general 
body  cases  with  six  motion  degrees  of  freedom,  there  are  21 
hydrodynamic  inertia  components. 6 These  are  best  given  in 
the  form  of  a symmetric  tensor.  Specifically, 

/j_  = If  2,  3,. ..6 

30i  0- ; dA 

an  : j = 1,  2,  3,...  6, 

where  0^  is  the  normalized  velocity  potential  caused  by  a 
unit  motion  from  some  reference  position,  and  it  can  be 
shown  that  the  various  d<t> i depend  only  on  body  shape. 

5n~ 

Miller^  has  shown  that  the  hydrodynamic  mass  values 
may  be  both  frequency  and  amplitude  dependent  for  irregular 
shaped  bodies.  For  this  condition,  the  set  of  differential 
equations  used  to  model  the  body  motion  become  nonauton- 
omous.  The  techniques  described  in  the  body  of  the  report 
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may  be  useful  in  obtaining  approximate  solutions  to  the 
types  of  equations  used  to  model  this  problem. 

The  effects  of  hydrodynamic  inertia  terms  on  some 
physical  systems  may  be  significant  and  should  be  considered 
when  writing  equations  of  motion  for  such  systems.  If  they 
are  determined  to  be  negligible  as  applied  to  a particular 
system,  they  should,  by  all  means,  be  neglected  in  order  to 
effect  solutions.  However,  in  the  case  of  certain  body 
motions  in  a viscous  fluid,  the  equations  of  motion  may 
become  highly  nonlinear  due  to  the  effects  of  damping  and 
hydrodynamic  inertia.  For  this  case,  numerical  methods  of 
solution  may  be  required. 
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